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ABSTRACT 


This thesis examines the combination of three algorithms: Graph Theoretic 
Tracker (GTT), Hough Transform, and Heuristic Search to enhance the detection of 
spectral tracks of underwater targets in LOFARGRAMS. Previous studies examined 
these algorithms separately. Here, GTT 1s used as a pre-processor of the LOFARGRAM 
display data to locate optimum paths of signals through noise. The line tonals in the 
output image from GTT are then manipulated by the Hough Transform into clusters of 
points in parameter space. A Heuristic Search sorting technique 1s employed to 
aepetmine cluster centers. These cluster centers are then reconstructed back into line 
tonals using the inverse Hough Transform formula. The results of this thesis show 
improvements by taking the Hough Transform of the original LOFARGRAM masked 
by the output of GTT. The effect of background noise ts offset by the accumulation tn 


the parameter space. Subsequently, the recovery of desired tonals 1s improved. 
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I. INTRODUCTION 


A. BACKGROUND 

As sonar systems continue to increase in size and 
complexity, operators are in danger of being overwhelmed by 
the information presented in displayed data. There is a need 
fem algorithms that “can “support automatic detection “and 
Mmmdaekinig Of Sigqnels of interest. In military applications, 
improvements in the area of detection of signals in noise are 
always of continuing interest. Concerning Anti-Submarine 
Warfare (ASW), one classical method of displaying passively 
received acoustical data is the LOFARGRAM. LOFAR is an acronym 
for Low Frequency Analysis and Recording. In the form of a 
waterfall display, acoustical data is presented spectrally 
with the y-axis representing time and the x-axis representing 
frequency. In the frequency domain, man-made signals possess 
a higher power content when they are compared to background 
noise. Because the shading of the pixels in the LOFARGRAM is 
power content dependent, pixels with higher power appear 


Beighter. 


B. PREVIOUS WORK 
A broad spectrum of techniques have been applied to 

detect acoustic signals passively in ocean noise. Two 
previously studied algorithms include the Graph Theoretic 
Tracker (GTT) and the Hough transform. Both algorithms have 
been examined in detail using synthetic test sets. Ross [Ref. 
1) gives a detailed analysis of GTT while Wang's study as 
devoted to the Hough transform (Ref. 2]. The main advantages 
of these algorithms include: requiring no @ priori knowledge 
of location or numbers of signals to detect; having the 
ability to detect single, multiple, crossing, and swept 
tonals; producing output suitable for immediate display; and 
providing accurate frequency estimates. The study here is 


aimed at improving the overall performance of tonal tracking. 


C. THE DESCRIPTION OF THE PROBLEM 

This thesis investigates the application of the GTT and 
Hough transform algorithms on a real data set provided by the 
Naval Research Laboratory (NRL). The use of synthetic data 
sets in previous studies allow tighter controle. 
experimentation and is therefore well suited for quantitative 
analysis. However, in this thesis, the performance Of chee 
algorithms using real data provides a more accurate 


qualitative evaluation. 


A tandem combination of GTT and Hough is used to maximize 
their respective advantages and minimize their inherent 
disadvantages. A final step involves the use of a heuristic 
search sorting technique. This step determines which clusters 
mimetne parameter space of the Hough transform should be 
reconstructed back into line tonals in the feature space. 

BEAMOO.BIN is a LOFARGRAM provided by NRL. Displayed in 
Figure 1, BEAMOO.BIN is a 256x256 byte image of three stable 
tonals and one swept tonal at a low signal-to-noise ratio 
EONR). In order to keep this thesis at the "UNCLASSIFIED" 
evel, both the origin and frequency content of this LOFARGRAM 


are not discussed and in fact are unknown to the author. 


7 nod aT we pra gies tay, 

Ty ‘so my SL LAE SZ 
ee aed oe Sine Pees 

re ne tae ces = aS 

r ae eek te coe ' een 


- aes i 
xu a a. Sac 





Theedirection of this thesis is best wl]llustraveq=b aa. 
flow graph presented in Figure 2. Because the GPP alooriemm 
Operates most efficiently with normalized data _ sets, 
BEAMOO.NRM, the normalized version of BEAMOO.BIN, is used as 
the input LOFARGRAM into GTT. 

The output file from GTT, TRACKS.B, is then ™%tilizedmae 
input for them@tough™ transtoun, After the image has been 
transformed into clusters in parameter space, a heuristic 
search sort is performed to determine which clusters are 
Suitable for reconstruction. Finally, the inverse of the 
Hough transform is employed to convert selected clusters back 
inco. line “ronda lsz 

This thesis consists of five chapters. Chapter I provides 
an introductory description of the track detection problem. 
Chapter II presents a description and implementation of the 
GTT algorithm. Chapter III emphasizes the Hough transform 
method. Chapter IV details both the original and improved 
heuristic search algorithm. Finally, Chapter V includes 


conclusions and recommendations 2c gesuat hemesroay. 
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Figure 2. Procedure Flow Graph in this Thesis 


Li. GRAPH THEORETIC TRACKER 


A. THEORY 


Graph Theoretic Tracker (GTT), developed by L.J. Wu and 
R.A. McConnell at Naval Research Laboratory (NRL), Washington, 
D.C. , is based on graph partition theory [{Ref. 3]. Jensen 
developed this optimum network partitioning theory whereby a 
set of elements are partitioned in such a way that a pre- 


defined cost function is optimized [{Ref. 4]. 


B. DEVELOPMENT OF THE ALGORITHM 
Applied to LOFARGRAM analysis, the track detection 

problem 1S simply translated into a graph partitioning 
problem. Optimum partitioning 1S utilized” co exyeraee 
features such as prominent line tonals. Each time line of 
the LOFARGRAM is mapped onto a graph where the individual 
pixels or frequency bins of the image correspond to nodes in 
the graph. According to Wu and Curtis [Ref. 3], 

...the partitioning is accomplished through the or@eaim 

enumeration of all possible partitions of agenaon 

followed by a recursive search using dynamic programming 

methods. The final partition generated by this algorithm 

is optimal with respect to some objective cost function 


used to drive the dynamic programmingescamen. 


Sample partitions of a graph are shown ine aagi-ce = 








—— a er 


eGie igre 3) 6 Meshes) ls hmalesa (ene 
LOFARGRAM to GTT 
Graph 


In order to maintain connectivity requirements, two dummy 
nodes are used as end nodes. Line tonals in the LOFARGRAM 
will appear as cuts through the graph. In order to detect 
these tonals, the optimum partitioning algorithm utilizes a 
mest function based on the signal-to-noise ratio (SNR). The 
memse EStimate needed for this cost function is found by 
Calculating the mean weight of the nodes between pairs of 


meaekS Or cuts. This is shown graphically in Figure 4. 





Figure 4. Cost hirer sen. iim. 


Symbolically, the cost function is defined as folle@wee 


CLiy = YOqy - Nye (278) 
where the Signal estimate, n., 1s determined by integqrauema 


the graph weights along the track path, t and the noise 


1/ 


estimate, oa.:. 


ie weighted by a scalar constant, y, is obtained 


by calculating the mean weight of the nodes between t= tii 
the scalar constant, y, can be used as a Means Gou—aee 
detection threshold, “Phe cose runerren, cfi;, 1s simply the 
cost of a particular pair of cuts through the @ocace By 


summing all possible cuts, the total cost can be determined: 


totalcost = )) cfy,. (2a 
i 





Figure 5. Various Possible Tracks 


C. IMPLEMENTATION OF GTT 
1. Parameters Involved 

Because the number of potential signals in a typical 
LOFARGRAM can be quite large, graph partitioning processing 
mmes Can become unreasonable. By limiting the number of time 
lines of data processed, GTT can produce output in an 
acceptable time period. This is achieved by limiting the 
height and width of the processing window. Currently the 
variables used to control this window include K, L, and P. K 
establishes the maximum number of pixels or frequency bins 
that a track can deviate from one time line of data to the 
next. A maximum value of four is allowed for K indicating 
that a track has nine possible positions in the succeeding 
mae, {-4,-3,-2,-1,0,+1,+2,+3,+4}, relative to its current 


position. 


10 


The parameter L determines the height of the 
processing window, limiting the number of lines of data to be 
partitioned. A maximum of three lines is permitted. E 
defines the width of the processing window and represents the 
number of frequency bins or pixels per time line. The maximum 
value allowed for this parameter is 512. By using these 
constraints, the number of possible cuts, N, evaluated reduces 


EO 


N = P(2K+1)4-1 ¢ <¢ QUP-2, (2535 


For comparison, Table 1 shows the number of cuts 
analyzed for a 2x256 window of data from BEAMOO.BIN with the 
number of cuts performed if the maximum values for all 


parameters is used. 


Table 1... ‘COMPARISON OF GUUS 


2X256 window of data maximum values allowed 
P = 256 P = 512 

Kee Ke 4 

L = 2 Le = 


xX 
il 
~ 
OV 
CO 
z 
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41472 


fags 


2. Source Code 
The software implementation of GTT, provided as 
Appendix A, contains the following programs, listed in the 


order in which they are used. 


@ MAIN.C 
@ PART.H 
@e INITIAL.C 
@ GENNODE.C 
See e DATE. C 
SreoARTIT.C 


@ LIBRARY.C 


In order to take full advantage of modularity inherent 
to the C programming language, the GTT algorithm is 
Baolemented aS a@ main program, MAIN.C, calling five separate 
source files. PART.H is provided as a header file to 
establish definitions of and limitations on the variables 
used. A graphical depiction of GTT activities is shown in the 


Plow chart provided in Figure 6. 


Ly 


oe —— 2. = ce 


Read in image ( initial.c ) 


Build tree ( gennode.c ) 


read LxP window of 
pixels into weight ( update.c ) 
array 


establish weighted | 
noise and signal ( partit.c ) 
estimate arrays 


partition graph 


output tracks.b (library.c ) 


Figure 6. GTT Flow Chart 


Using command line inputs, INITIAL.C determines if the 
parameter K, L, and P are within required limits and reads the 
LOFARGRAM in image format. GENNODE.C is then called to build 
m@me tree that will be graph partitioned. This module maps 
individual frequency bins in the display image onto nodes in 
eegraph. In the main loop of MAIN.C, GTT successively 
evaluates an LxP window of data as a graph to be partitioned. 
Within this window, GTT attempts to maximize the total signal 
while minimizing the noise. The output from this window is a 
set of cuts that represent the best potential tracks. Within 
Ene main loop, UPDATE.C is called to read this LxP window of 
pixels into a weight array. Graph weights for the nodes are 
determined by the pixel intensity representing those nodes. 
UPDATE.C then constructs arrays for the weighted noise 
estimate and signal estimate. When this 1S completed, 
meerltT.C is called to partition the graph based on the cost 
function (see equation 2.2). Finally, LIBRARY.C receives the 
Meee partitioned graph and outputs this as TRACKS.B. The main 
loop of MAIN.C is traversed again using the next LxP window of 
mata from the input file. Tiemowtmpweme tle,  TRACKS.B, is 
displayed in the same format as the input LOFARGRAM, which is 
a waterfall display of frequency versus time. Contrary to the 
input file which is displayed in 256 gray levels, TRACKS.B is 
binary with the black background representing the absence of 


Signal and white portraying the presence of signal. 
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For clarity in display, the inverse of the imagewes 
TRACKS.B is used in the following figures. The GTT output of 
BEAMOO.BIN, shown in Figure 7, was obtained with the parameter 
settings of K=1, L=2, and P=256. Although the GTT output of 
BEAMOO.BIN displays a large number of false alarms or noise, 
the three constant tonals and the swept tonal are visibly 


evident. 





aC 7, GTT Sco 
TRACKS.B, with 
y=6 
3. BEAMOO.NRM Image 
The LOFARGRAM BEAMOO.NRM is the normalized version of 
BEAMOO.BIN. BEAMOO.NRM is used as the input data set for GiTT? 
To determine the pattern of normalization, the IM-4000 Image 
Manager is used to analyze the histograms of both BEAMO0O0O.BIN 
and BEAMOO.NRM. Containing the full spectrum of gray scales 


0-255, BEAMOO.BIN possesses a mean pixel value of 161. 


iS 


BEAMOO.NRM however is normalized to 64 gray scales with a mean 
Bixel antensity of seven. in order to reproduce BEAM00.NRM, 


the following standard normalization equation is used: 


25 


Se ae 
Mee Ome mn pixel 


Norm Imag = Old_Imag~* ( 


Gir htewever Goes MOE, LUNCEIOn properly with the 
normalized image. Further study reveals that BEAM0O0O.NRM is 
probably constructed by normalizing each record or time line 
mepadta instead of global normalization of the entire image. 
Pmetead of analyzing this normalization further, this thesis 
concentrates on manipulating the output of GTT, using 
BEAMOO.NRM as input. 

4. Detection Threshold 

Mentioned earlier, the detection threshold can be 
memrpulated by changing y in the cost function calculation, 
equation (2.1). This can be achieved by numerically changing 
mee scalar constant used in the function, COSTFUNCTION(i,j), 
Of PARTIT.C. Figures 8, 9, 10 , and 11 show the effects of 


Manying y from five to ten. 
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Figure 8. Effect of Thresnela 
V2 





Figure 9. EIte€ce Creme esiorm 
Yee 


iy 
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Figure 11. Effect of Threshold 
y=10 
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As seen in these figures, the decision to set a 
detection threshold, y , at a certain level involves trading 
off detection for false alarms. If the threshold is set at a 
low value, as it is in Figure 8, the signal becomes more 
evident at the expense of introducing a large number of false 
alarms. As the threshold is increased to values of six, seven 
and finally ten in Figures 9, 10, and 11, the number of false 
alarms diminishes but the probability of miss increases. From 
visual observation, a detection threshold of six provides the 
Dest outpuc. 

The output of GTT shows the detection of prominent 
points along the track, but those points tend to be 
discontinuous. The next processing step involving the Hough 
transform is intended to detect potential line tonals existing 
in the output of GTT. The Hough transform achieves this by 
determining which disconnected points exhibit collinear 


tendencies. 


IPs, 


III. HOUGH TRANSFORM 


A. THEORY 

One method of detecting the presence of line and curve 
segments in images is the Hough transform. If a line can be 
mathematically defined by the equation y = mx + b, then the 
infinite set of points that comprise this line all possess the 
same value of slope m and intercept b. The Hough transform 
mses this fact to map collinear or nearly collinear points in 
the (x,y) plane of the image space or feature space to the 
(m,b) coordinate of parameter space. One obvious discrepancy 
occurs when the line assumes a vertical orientation resulting 
imi the slope approaching an infinite magnitude. This 
eenguilarity can be avoided by using the normal 
parameterization of a line first suggested by Duda and Hart 


(Ref. 6]. The equation then appears in the following form: 


9 = xcoso + ysine , Dour: (3.1) 


Using the parametric equation of a straight line, each 
collinear (x,y) coordinate will map to a sinusoid in parameter 
space. After transformation, an infinite number of sinusoids 
are generated, each intersecting at the point in parameter 


Space corresponding to the slope and intercept of the line. 
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1. Hough Transform Algorithm 
The Hough transform adheres to the algorithm displayed 


in’ Table 2 [Rema 7 ie 


Table 2. HOUGH TRANSFORM ALGORITHM 


For each (x,y) € P do 


0,7, a8 do 


For 6 
p = xcos6 + ysin®@é 
H(8,p) = H(0,p) + 1 

end do 


end do 


2. Properties and Illustration 
Using the normal parameterization of a line, the Houa@ 
transtorm method yields four main properties || Remger 
@e A point in the feature space corresponds to a sinusoidal 
curve in the parameter space. 


® A point in the parameter space corresponds to a linewema 
the feature space. 


@e Points lying along a line in the feature space correspond 


to sinusoidal curves through a common point in the 
parameter space. 
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Cro Miecu Vimcgmcl Onguac Simnusoidd! GCurvesin the parameter 
SeieemeColeeccooNeseo Lines throug a common point in the 
feature space. 


These properties are demonstrated using the illustration shown 


i Figures 12 and 13. 





Figure 12. Test Image for the 
Hough Transform 


ae 


ee, 
a 
<< 
—_ 
— 
= 

D> 
= 
= 

~~ 
———. 


>o© C&S 
theta Cdegrees) 





Figure 13. Sinusoids in Parameter Space from Collinear 
Points 


The line segment, displayed in Figure 12, is uUSéqmae 
an bese. Samage= The three indicated collinear points are 
extracted and used to generate the three sinusoids displayed 
in Figure 13. The intersection of the these three sinusoids 
in parameter space demonstrates the attractiveness of the 
Hough transform in locating collinear points in feature spaces 
One only has to determine those points in parameter space 
where a large number of intersections occur to locate lines in 
feature space. The parameter space iS viewed aS a two 
dimensional array constructed by quantizing p and @ Same 
cells. Each cell is assigned an accumulated value. Equation 
(3.1) maps each collinear point to a sinusoid 1n parameéwem 
space. Cells lying along the resultant  sinuse teas 


incremented by one. 


a5 


After the image has been transformed, accumulator 
aes HNavingmenugheecounts iandicate intersecting curves in 
Parameter Space and therefore lines in feature space [Ref. 7}. 

3. Reconstruction 

After transformation, the accumulator array is 
evaluated to locate those cells possessing high counts. The 
problem of reconstructing the cluster center (p,8) is solved 
eye performing the inverse of equation (3.1). Solving equation 


m1) for x yields 


sin68 


SOE ee) 





x=p - yl 


In order to implement this algorithm in a computer program, a 
Beterence point located at the center of the image is needed. 
mega Lelerence point simply introduces an offset (x,,y,) 


mero equation (3.1): 
p = (x-x,)cos0 + (y-y,) sind. (3-3) 


Solving this equation for x yields 


x=x a (354) 


B. IMPLEMENTATION 
The Hough transform is implemented using HOUGH.FOR, a 
FORTRAN program written to take advantage of the points 


mentioned in Section A.2. 
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Previous programs of the Hough transform in [Ref. 8] and 
in [Ref. 2] involve the detection of signals present in noiae 
and are therefore gray scale oriented. Because the output of 
GTT is binary, a new program, HOUGH.FOR, is used which is 
provided in Appendix B. The flow chart in Figure 14 depicts 
the flow of operation from output of GTT to tama 


TECONStruction. 


read in ( gttconv.for ) 
TRACKS.B 


establish 
reference point 
Xo0,Yo 


Hough 
Transform 


increment 
accumulator 


threshold 
accumulator 





Figure 14. Hough Transform 
Plow Ghar 
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1. Source Code 
Becausesene output of GTT is a fixed length 128x512 
fasertile, a pre-filter 1s needed to convert TRACKS.B into a 
Meoxz256 byte array. Ine@rORDRAN program, GITCGONV.FOR, that 


performs this operation is listed in Table 3. 


Table 3. GTTCONV.FOR LISTING 

ee 
program goutconv 

THis pregram takes a fixed Length 128x512 byte 
output file from the GTT algorithm, which is 

in binary format, and converts it into a 256x256 
byte array. The output file, OUTFILE.DAT, is 

now suitable for input into the Hough Transform 
eaeigorithm: 

Pyiecuol image, 120,512) 72> image 256,256) 

imMeceder a image, 256,256 ) 
open(unit=1,name='tracks.b',status='old', 


eo Goo 


ze AGCGecS=—Glmecece ~reCconds1ze=l2e8ymaxrec=128) 
open(unit=2,file='outfile.dat',status='new', 
* access='direct',recordsize=64,maxrec=256) 
@o 1=1,128 

BeAC(Marijiol IMagqe(i,9),  j=1,512) 

de je 5.6 

b2_ image((i*2)-1,})=bl_image(i,j) 
end do 


do qW=257e5 1 2 
b2 image({1*2,)-256) pl image 1, 4) 
end do 
end do 


do i=) 256 
Write(2ame) (D2 lmageeignr a), 250) 

end do 

close(unit=1) 

close(unit=2) 

end 


Zo 


Subroutine h in HOUGH.FOR performs the Hough 
transform. 8 is incremented from 0 to m radians in intervals 
of 1/256 radians. For each nonzero pixel in the image, ene 
Hough transform equation is executed and each cell in the 
accumulator array that is traversed by the resulting sinusomd 
is incremented by one. After performing the Hough transform, 
the accumulator array 1s thresholded and normalized for 
display on the IM-4000 Image Manager image processing system. 
This image processor requires that the scale for normalization 
range be from 0 to 243. Levels 244 to 255 are internally 
reserved (Ref. 9]. Cells in the accumulator array greater 
than or equal to a user defined threshold will map to the most 
likely set of collinear points in feature space. Subroutine 
R in HOUGH.FOR is then called to reconstruct lines from 
cluster points in the accumulator array. After determining 
the number of cluster points and their location, (p,8), the 
reconstruction equation (344) tis wsed. 

2. Output 
a. Test.dat 
The test image, TEST.DAT, shown in Figure lS) 
used as the input file for the Hough transform. The result Gf 


this transform is shown in Figure 16. 
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Figure 15. Test Image in 
Space Domain, 
VES. DAT 





Pyroguec 6 Hough Transtorm Of 
eo OAT 
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As predicted, the sinusoids intersect at a common 
(p,8) location. The reconstruction of this Hough =transtemaee 
image is displayed in Figure 17, which demonstrates that in 
the absence of noise the Hough transform can faithfully and 
easily reproduce the line test image. The accumulator array 
for the Hough transform contains one cell that possesses the 
highest number of sinusoids traversing it, and therefore the 
highest count. In three dimensions, this appears as a lone 
peak rising from a plane defined by p and 8. The 3-D plot of 


the Hough transform in Figure 16 is shown in Figure 18. 
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intensity 


Figure 


ee 





Figure 17. Reconstructed Image 
of TEST]. DAT 





S-OeelOt Of Ehe Hough Transform 
of TEST.DAT 
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b. fTracks.b (Outfile.dat) 

The output of GTT, TRACKS.B, shown in) Pigure yee 
converted, as discussed previously, into OUTFILE.DAT anda. 
as input for HOUGH.FOR. The Hough transform of this file is 
displayed in Figure 19. Because TRACKS.B contains numerous 
false alarms in the form of small line segments, the @hog@ee 
transform appears to be extremely noisy and convoluted. After 
the accumulator array for this transform is thresholded at a 
value of 100%, reconstruction produces the singlé most 


dominant tonal in the LOFARGRAM, which is shown in Figure 20. 





FLQULES™ 1 Senin G bie Lilet Serr lt mmeohts 
GTT Results, 
TRACKS.B 


SL 


Figure 20. Reconstructed 
Dominant Tonal, 
with 100% 
Threshold, of 
Figure 19 


The reconstruction subroutine in HOUGH.FOR did not 
sense the less evident tonals because the accumulator cells 
associated with those lines contain counts much lower than the 
@ell linked to the most prevalent tonal. This is evident in 


Figure 21 which shows the 3-D plot of Figure 19. 
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Figure 21. 3-D Plot of the Hough Transform 
of GTT Results, TRACKS.B 


Even with a threshold value of 70%, the most 
dominant tonal is still the only one reconstructed, as shown 
in Figure 22. Clearly then, an alternative is required to 
extract those cluster centers associated with tonals of 


interest in the LOFARGRAM. 
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Hagure 22. Reconstructed 
Dominant Tonal, 
with 70% Threshold, 
of Figure 19 


The issue is how to find relevant cluster centers 
in the parameter space of the Hough transform. In the next 
chapter, a heuristic (greedy) type of algorithm is developed 


to locate cluster centers. 
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IV. HEURISTIC SEARCH 


A. BACKGROUND 

Clearly a simple thresholding operation is insufficient to 
reconstruct all applicable tonals from the accumulator array. 
A method needs to be developed to effectively sort the 
accumulator array and reconstruct lines from accumulator cells 
of interest. Two methods, LAS cluster technique and sorting 
technique, have been previously studied as a means to perform 
cluster analysis [Ref. 2]. The Land Analysis System (LAS) is 
a software package encompassing a variety of functions 
designed to process and analyze image data. Based on a 
Digital Equipment Corporation (DEC) VAX 11/780 computer, the 
LAS runs under the VMS operating system | ker ace Three 
programs of interest within the multispectral processing 
functions of LAS include HINDU, ISOCLASS, and KMEANG The 
HINDU program performs an unsupervised classification based on 
multidimensional histograms. ISOCLASS performs an unsupervised 
cluster classification and KMEAN groups input image pixel 
values into a predetermined number of clusters. In addition 
to applying a sort based on threshold, Wang studied these 
three programs as an alternative cluster analysis method. 
This thesis examines a possible third procedure, heurisiae 


search. 
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B. THEORY 

The heuristic search method involves scanning the 
e@eccumulator array and determining the location of groups or 
Clusters of cells or points which possess values equal to or 
greater than a user defined threshold. The threshold is based 
on percentage with 100% assigned to the cell with the highest 
@ount. As the search proceeds, the clusters are condensed 
mimeo a specific eGluster or cells which will then be 
reconstructed. In order to implement this procedure, a cost 
m@onction is defined. 

1. Cost Function 

The LOba! Cem rimecer on 1s teommed from the addition of 

mecost. function based on the number of points (cost,) and a 
cost function based on distances between points (cost,). 


Symbolically, these functions are illustrated below. 


cost, = K,N (4.1) 

gost, -"KD= ) >) ku, — x, (4.2) 
Uy, Xy 

totalcost = cost, + cost, (2238) 


mgere K, and K, are scalar constants used to weight cost, and 
cost,, respectively, u, represents cluster group i and x; 
Meeeesents the fj points within cluster group i. An 


mebustration shown in Figures 23 through 26 is used to 


demonstrate the cost function. 


36 


Figure 23 shows an example of an accumulator array 
with cells containing values at or above the threshold 
depicted as black points on the grid. Numbers’ listed 
vertically mark rows in the array; horizontally listed numbers 


mark columns. The black x's indicate cluster groups. 





Figure 23. Test Image of 
Heuristic Search 
With N=4 wo=O0e0 


At this point in the calculations, N equals four and 
a marker 1s placed on each point yielding a distance from each 
marker to its respective point zero. The minimum distance 
between markers 1S now computed, and a new marker is placed at 
the mid-point between the two closest markers. In Figure 24, 
N now equals three, and the distance cost associated with the 
new marker is calculated. The distance used in this 
calculation is the distance from the new marker to each of the 
Original cluster markers. At the completion of each Stdtey 


the total cost is computed and tabulated. 
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Figure 24. Test Image of 
Heuristic Search 
with N=3 D=5.0 
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Figure 25. Test Image of 
Heuristic Search 
with N=2 D=11.32 
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Figure 26. Test Image of 
Heuristic Search 
with N=1 D=34.81 


When all states have been evaluated, the state 
possessing the minimum overall cost is reconstructed back into 
line tonals in feature space. The three cost functions are 
shown in Figures 27, 28, and 29. The total cost function for 
the previous example is plotted against states one through 
four, 

It is evident that state two representing the two 
cluster group situation is the lowest cost configuration which 


would be the answer for the number of reconstructed clusters, 


BY 
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After experimentation, scalar constants K =1.0 Jana 
K ,=0.1 were chosen to ensure that, within the total cost 


function, cost, 1S properly eet es: ely eecosn.. 


Cc: IMPLEMENTATION 
1. Source Code 
The Heuristic Search method, written in the FORTRAN 
programming language, is provided in Appendix C as HS.FOR. 
This program takes the 256x256 accumulator array constriu@ged 
in HOUGH.FOR as input and determines the minimum cost 
configuration. This configuration 1s then provided as igus 
for RECONST.FOR which reconstructs the optimal minimum cost 
configured accumulator array back into lines in feature space. 
a. HS.FOR Program 
It is important to remember that the image 
supplied to HS.FOR has already been thresholded by HOUGH.FOR. 
Pixels containing values equal to greater than the threshold 
are retained and those below the threshold are set to zero 
(pixel value 0 equates to the color black for display 
purposes). 
After conversion from a 256x256 byte image into an 
integer format, the x,y locations of the N non-zero pixels are 
stored in four 1xN arrays: S PEAK ROW, S PEAK COL, © PEAKRRe@se 


and C_ PEAK COL. 
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S PEAK _ROW and S_PEAK COL store values from the 
input sample image and remain unchanged as the reference 
image. € PEAK ROW and C PEAK C@L represent the cluster image 
and store the locations of the current cluster groups. The 
length of c_ peak_row and c_peak_col arrays will change as the 
program progresses from state n=N non-zero pixel or peaks to 
State n=l. Initially, these four arrays are identical. 

Within the main program loop, subroutine CENTROID 
performs the majority of the processing on the cluster image. 
Sequentially, this subroutine determines the minimum distance 
between peaks and places a new marker at the mid-point of the 
two peaks involved in the minimum distance calculation. It 
adjusts c_ peak _row and c_peak_ col arrays to include the newly 
calculated midpoint and determines the distance factor to be 
used in the distance cost portion of the total cost function. 
It then performs the next set of distance calculations as the 
loop is executed again. After all states have been evaluated, 
Subroutine min_cost is called to determine and output the 
cluster group configuration representing the minimum total 
eOsSt. 

b. RECONST.FOR Program 

This program takes the output from HS.FOR as input 

and performs reconstruction calculations nearly identical to 


those found in HOUGH.FOR. 
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The main difference is that the pixels )Seaume 
reconstructed now represent the cluster configuration Ome 
minimum cost. 

2. Testing 
a. Threshold 30% 
At a threshold of 30%, ten cells dinars 


accumulator array are encountered and are listed in Table 4. 


Table 4. ACCUMULATOR ARRAY CELLS OF 
HOUGH TRANSFORM PARAMETER Vor neo 


ee 


254 ee 6 
250 Z6, 


ie 
| 8 | 129 
—— a 





After reconstruction from parameter Spaceniae 
feature space, these ten peaks yield ten overlapping lines 
shown in Figure 30, depicting the most dominant tonal (iis 


Original LOFARGRAM. 
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Digumeme0 -  Recomstmucted 
Dominant Tonal with 
30% Threshold of 
Figure 19 


In order to implement the heuristic search 
Cluster analysis, the image shown in Figure 30 is used as the 
input image for HS.FOR. The ten lines correspond to ten 
States initially. As the program proceeds from state n=10 to 
eeate n=l, cost,, cost,, and total cost are computed. Cost 


function values for this example are listed in Table 5. 
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Table 5. COST FUNCTION VALUES FROM 
FIGURES 27m 28, AND 29 


state, cost, cost, total 
cost 
30.0 | 0 0.0000 nO oad 
a 
a 
ae 
ae 
ee 








6.8284 4.6828 





8.2426 300 es 










10.2426 50243 










I26> 058 1270S 









It is evident that state n=2 represents the 
minimum cost configuration. Two peaks in the accumulatom 
array result in two lines in feature space as shown in 


Paguine so 4M 
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Figure 31. Reconstruction 
of Figure 30 
after Heuristic 
search 


b. Threshold 20% 
When the threshold is reduced to 20%, the number 
of accumulator array cells or peaks meeting the threshold 
Criteria iS increased to 32. The 32 lines reconstructed 


from these pixels without heuristic search are shown in 


maigure 32. 





Higuse 32 Reconstructed 
Tonals with 
20% Threshold 
Gnerigune 19 
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The sloping line to the far left in chews 
represents the top portion of the slowly varying curved 
tonal in the original LOFARGRAM. The next set of vertical 
lines depicts the less evident tonals while the majoricyees 
lines are centered around the position of the domunmame 
EOnal The two vertical lines to the far righty 
Spurious and represent noise pixels that met the threshold 
criteria in the accumularonvarray. 

The heuristic search yields a minimum cost state 
of eight cells, and the reconstruction of which is Showa 
Figure 33. The large number of lines near the dominant 
tonal has been reduced to three. It is interesting to note 
that the two spurious lines associated with noise are 


retained. 





Figure 33. Reconstruction of 
Bigune (S2anberm 
Heuristic Search 
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c. Threshold 19% 
Lowering the threshold to 19% yields 40 peaks in 
the accumulator array corresponding to the 40 lines shown 


in the image in Figure 34. 





Figure 34. Reconstructed 
Tonals with 19% 
Threshold of 
Figure 19 


New additions from the 20% threshold image include 
Miemvertical Line to "enemitar Went representanmg another 
Spurious noise cell in the accumulator array. It also 
includes two more sloping lines depicting the mid and lower 


portions of the swept tonal from BEAMOO.BIN. 
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Eleven cells represent the minimum cost state in 
the accumulator array resulting in 11 lines reconstru@ees 


in feature space as shown in Figure 35. 





Figure 25. Reconseruceron on 
Figure 34 after 
HeuUpiStt Gaecaren 


The spurious noise lines have unfortunately been 
preserved while the majority of lines clustered around 
the dominant tonal position have been reduced to three. 

It is evident that the heuristic search cluster analysis, 


although working well, needs to be improved. 


49 


3. Improvements 

In order to reduce the effects of noise while 
Betaining tnose cells @am the aqecumulator array associated 
with the signal tonals, the input image for the Hough 
EaansrOnm 1S Subjected fe further processing. MULTIPLY.FOR, 
provided in Appendix D, is used to multiply the output of 
GTT, TRACKS.B, with the original LOFARGRAM, BEAMOO.BIN. The 
resultant image is a replica of TRACKS.B with the 
exeception that every non-zero pixel contains the value of 
BEAMOO.BIN at that (i,j) location.This image now emphasizes 
imeesional tonalss Phe output of MULTIPLY.FOR, OUTIMG.DAT, 
is now used as the input image for the Hough transform. 
Because HOUGH.FOR is binary image oriented, changes are 
made in the source code of this program to accomodate the 
O—255 Gray scale orientation of OUTIMG.DAT. These changes 
are listed as HOUGH G.FOR in Appendix E. 

a. Threshold 17% 

At a threshold of 17%, 38 peaks in the accumulator 
array meet the threshold criteria. When this array is fed 
into the heuristic search algorithm, HS.FOR, nine peaks are 
imetainead £Or reconstrucction.meine output of HOUGH _G.FOR is 
shown in Figure 36; the output of HS.FOR is depicted in 


Bargure 37. 
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Figure 36. Reconstructed 
Tonals with 17% 
Threshold of 
Figure 19 





Figure 372 RecOmseaverm onses 
Figuse 3654 e=cm 
Heuristic Search 
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Figure 37 shows that the manipulation of the input 
image by the multiplication of TRACKS.B and BEAMOO.BIN 
results in an output image from heuristic search where the 
main, secondary, and portions of the swept tonal are 
reconstructed without the introduction of noise. The 
threshold is then reduced to 16% resulting in 52 peaks 
meeting the threshold criteria. After the heuristic search 
is completed, 11 peaks are retained which reconstruct 11 
lines in feature space. These images are shown in Figures 


38 and 39. 





Figure 38. Reconstructed 
Tonals with 16% 
Threshold of 
Figure 19 


BZ 





RECONSEruCct won. on 
Figure 38 after 
Heuristic Search 


Figure 39% 


Again the main, secondary, and swept tonals are 
retained. However, the two vertical lines to the far right 
in Figure 39 are spurious resulting from two noise cells in 
the accumulator array passing the threshold criteria. 
With the proposed manipulation, it is obvious that more 
line fragments of desired tonals can be extracted out of 


Che -original image. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


A. GENERAL 

This thesis has presented an examination of the 
combination of three algorithms: GTT, Hough Transform, and 
Heuristic Search to enhance the detection of spectral tracks 
of underwater objects. Because signals of interest are best 
visualized using LOFARGRAMs, LOFARGRAM analysis remains an 
mueedral part Of tracking in the frequency domain. The 
combined approach of the application of these algorithms is 
used to take feel advantage of their respective 
characteristics. Conclusions drawn from this research are 
listed below for each of the three algorithms. 

1. GTT Algorithm 

@ Because the output is already in the form of potential 
tracks, it can be used for further processing. 


e This algorithm is able to distinguish between connected 
Signals and relatively poorly connected noise. 


@e GTT can process multiple stable tonals clearly. 
2. Hough Transform Algorithm 
@e This transform can extract desired line fragments of 


stable and sweeping tonals from planar point sets. 


® The Hough transform can handle spectral tracks buried in 
heavy background noise. 
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MuicealGotrrthm can Operate efficiently on both binary 
and gray scale images. 


3. Heuristic Search Algorithm 


@® Heuristic Search provides an easily understood 
alternative to LAS routines and sorting technique for 
cluster analysis. 


e This algorithm provides a search foundation for further 
image processing techniques such as simulated annealing. 


B. RECOMMENDATIONS 

Several unresolved problems remain at the conclusion of 
mais thesis. The cost function used in determining which 
cluster(s) to reconstruct in the heuristic search algorithm 
needs refinement. Presently, only the number of peaks and 
distances between peaks are considered. One) possi. Le 
improvement might include the addition of a height function 
which would favor those cells in the accumulator array with 
Mmrener counts. This could offset the effects of noise and 
emphasize the desired tonals in the LOFARGRAM. A second 
recommendation involves expanding heuristic search into the 
Simulated annealing algorithm. Simulated annealing is very 
eueeective in determining the global optimal solution and could 


be used in cluster analysis by finding the optimal minimum 


Bemmtion to the cost function. 


DS 


To operate effectively, the simulated annealing algcmeuaas 
needs to be given an initial guess of the location of 
potential tracks. The fragmented output from GTT is Not 
Suitable as an input state for simulated annealing. By using 
the Hough transform to convert the discontinuous powneemeaa 
GTT's output to line tonals, an initial guess can be provageg 
for simulated annealing which is close to the actual position 
of the tonals in the original LOFARGRAM. Research has bee@ 
conducted to apply the technique of simulated annealing to 
sonar track detection, but it only involved the usemiga 


laboratory generated synthetic data sets [{Ref. 11}. 
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APPENDIX A: GRAPH THEORETIC TRACKER SOURCE CODE 


MAIN.C 


#define MAIN 
#include <stdio.h> 
#include <time.h> 
#include <signal.h> 
genciude “part.h" 


Matn{( argc, argv ) 
omit argc; 
eoar *argv( |; 


{ 


Oe line=0, i; 

long starttime, stoptime; 

FILE ATitcialise(),. *anfile; 

VERTEX lastnode, GenerateNodes(); 

void Commoner ( ),eoaneition(), update); 


DumpNodes(), DumpTable(); 


[ERK KKEKKKKKKK KEK KKK / 


/* Boa@y of Program */ 

ee 

Gord) peinrre (ss \n\n  , verid); 

infile = initialise( argc, argv ); 
(ord Phiten bUdTGIng tree ...\n"); 
lastnode = GenerateNodes(); 

(void) printf("Nodes=%d\n\n", lastnode); 


i) 
line++; 
update( lastnode, infile ); 
if(feof(infile)) break; 


Veteimprintr ('Precessing group %3d”", line); 
(void) time(&starttime); 
Babe thloOn@lasenode): 
(void) time(&stoptime); 
stoptime ~-= starttime; 
(void) 
Beamer ("(%$21a:%21d)\n",stoptime/60,stoptime % 60); 
writepart( lastnode-l ); 
j 


i 


PART .H 


#define verid "ee * Pia tig. Vo aa 

#define MAXNODE 2500 /* max. tree nodes */ 
#define MAXPIXELS 256 /* pixels per linewey 
#define MAXLINES 3 /* number of lines */ 
#define nullset ~] 

#define TRUE 1 

#define FALSE 0 

#define BIGINT 10000 

typedef char BOOLEAN; 

typedef short VERDE, /* node in tree */ 

/* typedef int void; A /* dummy void type sim 


compiler lacking one */ 


/* Define Global Variables and Arrays */ 
#ifdef MAIN 


Shere int setsize[MAXNODE][{MAXLINES]; 

Shome wat setsizel{MAXNODE] [{MAXLINES]; 

unsigned int Ssize{MAXNODE];/* tree node sizes */ 
unsigned int acc(MAXNODE)];/* cost funct track dacoug 
ie table[{MAXNODE]; 

VERTEX opt [MAXNODE]; 

ae KLimit, lines, npix;/* Case uzen~ 
#else 

extern short int setsize{ ){MAXLINES]; 

extern short int setsizel{]{MAXLINES]; 

extern unsigned int size j7/7*=snedevsizc> a 

extern unsigned int acc[];/* c.f. track aces, 

extern int table[ ]; 


extern VERTEX CPE lg 


extern int KLimit, lines, npiusc 
#endif 
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INITIAL.C 


#include <stdio.h> 
menciude <string.h> 
mencilude "part.h" 


/* Parse Input Arguments */ 
POLE *initialise ( argc, argv ) 
ic Argc; 

ear *~argv{ |; 

{ 


enar imagefile[{[64], *strcpy(); 
FILE *infile; 
void Exe jeeisage(), Valermay, err()? 


iar gc > 1) 
Pircememe(aragumidy), /i: j==0;;argc>5) usage(argv(0)); 


Swarren( arge ) 


{ 


case 5: KLimit = atoi( argv[4] ); 

case 4 lines = atoi ( argv{3] ); 

case 3: Mpiee= waeOl( argy| 2] ); 

case 2: (void) strcpy( imagefile, argv[1] ); 
default: ; 


} 
Switch( arge ) 


{ 
CacSecmieww Vold) i plincnh( stderr,” input 
filename:\t\t" ); 
(vold) scanf("%s", imagefile); 
Case 2: (void) fprintf(stderr, “Line length:\t\t" 
)i 
Oy elcieSeant4 = sa. - ) Snp1x ); 
ccm VOTGnmeEDirIntnh(Staerr, Seqmentation 
imees: \t"); 
(void) scanf("%d", &lines); 
eicem- | VOld) tprintr(stderr,"Enter KLimit 
value:\t"); 
(void) Vseanfi("sd", &KLimit) ; 
} 


a9 


infile = fopen(imagefile,"r+b"); 
if (infile==NULL) err (“Unable to open input 
Pie 5 ye 


if (npix<l |; npix>MAXPIXELS) 
err("Invalid number of pixels"); 


if (lines>MAXLINES ,; lines<1) 
err("Invalid number of lines"); 


if (KLimit<] |) KE 4. 
err( “Invalid Kian, 


(void) Unlink( “trackoes ). 
Peturn( ines te). 


} 


VO1d uSsSage( argv ) 
Char “*argyv, 


(void) printf("\tusage: %s [input] [pixels] [{lines] 
Air yn aren, 
} 
void err (arg) 
Char arg, 


(VOId)y print t (45) 007 nace, 
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GENNODE.C 


#include "part.h" 
ert. node=0; 


VERTEX GenerateNodes () 


{ 
Begiseer Int 1, J; 
void Gen(); 
ye Iimitialise setsize */ 
mom ( 1=0; a CMAXNODE ; i++) 
fome( 1-0; «lanes. 74+ ) 
Setsizel(/ijij] = setsize([ij[j] = nullset; 
/* Generate tracks */ 
aoa (a — F< ee) 
{ 
setsize[node][0] = i; 
Gen(0O, 1); 
} 
/* Fill in node tables */ 
eOG (t=O 5 hayes: 1++) 
MOT wa) conoce..4 ++ ) 
if( setsize[j][{i] == nullset) 
setsize[j][i] = setsize[j-1l][1]; 
for (1=0;i<lines;i++) 
Poe (ai 0) MOCe. ++ ) 
setsizel[{j])[il]=setsize[j])[1i]+1; 
momen 1=0; 1<node: i++) 
size[i]}=0; 
mor ( t=O; 1<node; i+) 
for (j3=0;j<lines; j++) 
Size[i] = setsize[iJ][j] + size[i]; 
imei nOde ) } 
} 
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void Gen (LineNumber,BinNumber ) 
int LineNumber, BinNumber; 
{ 

register int ae 

Vou err(); 


if (node >= MAXNODE) err ("Too many nodes!"); 


LinéeNumbers+; 
for (1 = “KLimit 71<—hbere, 2.) 
{ 

] -= BinNunberara, 


Vf) > OS Rea <1 inate) 
{ 


setsize[node][{[LineNumber] = ]j; 

if (LineNumber == lines-l) 
node++; 

else 


Gen(LineNumber, }j); 
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UPDATE.C 


fac bua@e <stdio.h> 
meenclude “part.h” 


void update(lastnode,infile) 
int lastnode; 
PILE *infile; 


i 


register int ie eae; 
unsigned char buffer[(MAXPIXELS ]; 
Danie weight [MAXLINES]} [(MAXPIXELS}; 


Pom  j=-O0, Jcoheness) ate. ) 


if(fread(buffer,1,npix,infile)==0) return; 
Comet Cen pie; 1t+ 9) 
weight{j])[1] = (int) buffer[i}; 
} 


/* set up Signal estimate acc */ 
for (1=0;1<lastnode: i++) 
acc[{i]=0; 
for (1=0;1i<lastnode; i++) 
EOme( j=); fcrimes g++) 
acc[1]+=weight[j][{setsize[{i}[(j]]; 


/* set up table */ 
for( j=0; j<lines; j++ ) 
Om a en Ddix? 1++. °) 
weight (j](i]+=weight(j][i-1]; 


for( 1=0; i<lastnode; i++ ) 
maple | a j=1; 
for (1=0;1<lastnode; i++) 
Eom 7 Os <tiness y+ ) 
table[i]+=weight[j})[setsize[i][j]-1]; 
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PARTIT.C 


#inelude “pare. nh 

#include <stdio.h> 

VOlG Partitaon( nedeeoumes) 
int nodecount- 


{ 
EeGicter ant Loy} eee 
BOOLEAN Status; 
AGT cost[MAXNODE]; 
FILE *fp,*ftopen(). 
int mMinval; GCcoseringEerona., 
tor (1=0; 1<nedecount ra. ) 
Cose| 2 | = BIGING. 
Opti 1 || C= mt pice, 
fOr (1=17i1\medecoune, 7+) 
for( j=0; j<i; j++ ) 
Status=TRUE; 
forl(e=0- GClanes er. ) 
if(setsize[i][c]<=setsizel[j][{c]) 
Status-F ASE. 
break; 
} 
if (status) 
{ 
c = costfunction( i, je) sveeeeny | 
1LE( (EOSEl1) ) 7 - ee 
COS 1) Vasc; 
opt[i] = j; 
} 
} 
} 
} 
Ine EOS EUG EOn( @ ema) /* cost of seqnenen 
ERED 2 eee oe 
{ 


FEGiSter Int 1coce- 
Ane GOS 


icost = ¢int) (tablef[i]-table[j])*6; /*thresnoiee 
icost /= (int) (sizela |j=sizemioe 
reCurn(icost =" (ant) vacewaee 
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LIBRARY .C 


munclude <stdio.h> 
menclude “part.h” 


void writepart( node ) juries partitlon */ 
VERTEX node; 


{ 


register int 1, j; 
unsigned int bitmap[MAXLINES][MAXPIXELS]; 
PIL “Ouerale. 


Bom I Lanes 144 |) 
orem sO; me Ec, ee) 
Pileman | aia) = 0; /* bitmap O is 
gray scale black */ 


node = opt[node]; 
while( node != nullset ) 


mon ( i — Or biness i++ ) 

/* pixel quantization for the output 
image is 8 bits/pixel, with the 
gray scale running from 0 (black) 
to 255 (white) */ 

bitmap[1][setsize[node][{1i]}] = 255; 

node = opt[node]; 

} 
eerile = fopen(”tracks.b","a+b"); 
fOr (1=0;1<lines;i++) 
(Vote) eanialieotOLlemap| 1 ),l,npix,outtile); 


pvomenereclose(oucrile) ; 
REEVE; 
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APPENDIX B: HOUGH TRANSFORM SOURCE CODE 
HOUGH.FOR 


program hough 

byte b_ image(256),r_image(256) 

integer counter,u(256) ;max,* sine 567. 

integer gray._level, hl theta(2Z56 7 ieee cy 

integer i image( 256,256) , accum(256, 256), 1a eee 

real x0,y0,rho(256) ,*(256,256)) 917 See ee 

real theta(256),delta_ theta,rho max,rho 0,delta rho, facruam 


open(unit=1,name='test.dat',status='old',access='direct', 
*recordsize=64) 

open(unit=2,name='testh.dat',status='new',access='direct', 
*recordsize=64) 

open(unit=3,name='testr.dat',status='new',access='direct', 
*recordsize=64) 


Cc ***Begin Main Program*** 
Ee Read in test image data. 
Go 7=17256 
read(1'j) b image 
do \2=17,256 
1 image(i,j])=b_image(i) 
enddo 
enddo 
GOs =152 56 
dO. 1 =3F 7 5.6 
i1f(i_image(i,j).LT.0) i _image(1,})=1_1mage(i7 a) ee 
enddo 
enddo 
Ee FF NOUGI: Uieams © Orin 


call h(i_image,accum,theta,delta rho,rho 07.0 a 
Max-accum( 17 |) 
min=accum(1,1) 
Go. = 19256 
do, 1=1,256 
if(accum(i,})).«GT.max) max=accumenag. 
if(accum(i, }).LT.min) min=accun@ aap 
enddo 
enddo 
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Normalize accumulator array for display 
on IM-4000 Image Manager. 
factor=243.0/(max-min) 
do j=1,256 
Go 1=1 7572 56 
AeGumien | = Ninh accum( 1, )])-min)*factor ) 
Peccaeeunm i, | seGilet2 /) then 
b image(i)=accum(i,j)-256 
else 
b_image(i)=accum(i,j) 
endif 
enddo 
write(2'}) b_image 
enddo 


Me ECCONSERUGCEREDOM NOuUgh ECransSiorm* ** 
cma cacclil mc iaeNctasiernio,  one,delta rho, 


* ea @ mi OGe a0 ay OMe) C) 
doug — i, 250 
G@ = 12 56 


Pa xen me). 127) x int (i, 7)=x int(i,7])-256 
GaeiMaGe(ey— x Mt ( 1,7] ) 
enddo 
write(3'j) r_image 
enddo 
end 
‘see Main Program*** 


Lite RK AK A AKA KKK AAR KR KKK KKK KEKE KEKKEKKEKKEKKKEKKAKE 


subroutine h(i_image,accum,theta,delta rho,rho 0,x0,y0) 
dimension 1 image(256,256) 

integer accum(256,256) 

integer Gray level 

real theta(256) 

real x0,y0 

real pi/3.1415926/ 


Miata 2eaknhe accumulator array to zero. 
do j3=1,256 
Gomaa=—) 2516 
accum(i,j)=0 
enddo 
enddo 


DiCsecVletemeretaerrom 0 LO pl radians. 

delta theta=pi/float(256) 

do 1=1,256 
theta(1i)=float(i-1)*delta theta 

enddo 
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Q 


Fix the center of image@as srnemonta ane 

rho max=sqrt(float(256*256) 48 loat(256-25a0m 
rho_0=rho_max/2.0 

delta rho=rho _max7itloear(Zaqy 

“X0=fleac( 25677) 

yO=Lloat( 250/772) 


GO sige 6 
vaay, 
d6"ax=1 7256 
gray_level=i_image(ix,iy) 
if(gray level, LE. 0 \sqolvoms 
x=1x 
do: 1t=1)256 
Hough Transformation Equation. 
rho=(x-x0)*cos(theta(it))+(y0-y) *sin( theta (gee 
r=(rho+rho 0)/delta_rho 
ce jee) 
accum( it, ir) =accum(te, ue). 


enddo 
continue 
enddo 
enddo 
return 


end 


wake Kkeke ake KK KK KK KKK Kh KKK KKK KKK KKK KKK KKK KKK KKK KKK KK KK KKK 


Subroutine r(accum,max,h theta,h rho, rho, cel lama 
ENO 60; 0870 xe ee) 

integer accum(256,256),max,h theta(256), hi sne( Zany 

integer x _1nb (256,256) un 256) Goumecs 

real rho( 256), deltas che, aioe ewer. 70 

real. x(256,296)5 p17 3.14 lao 6, 


counter=0 
GOn je s.2 510 
do <= e256 
if(accum(i,j}).GE.max) then 
COUNEEr—Counter pe 
Determine the theta and rho locations 
in the accumulator array with the 
highest count value. 
h_ theta(counter)=i 
h_ rho(counter)=}j 
endif 
enddo 
enddo 
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do 3=1,counter 
GiO@gy—mmeeno( )})*delta rho-rho 0 
enddo 


Gow j—1, Counter 
GOr =e 56 
Reconstruction Equation. 
x(i,j})=x0+ 
(Ne (p(y 0-1) *Sami (i theta(;)—1)*pi/256))/ 
COs Ny Eneta( ) )- iyap17 256 ) 
x_int(i,j)=jnint(x(i,j)) 
enddo 
enddo 


Go” j=1,256 
do l=1,counter 
Tel) cee mena), 1) 
enddo 


do m=l1,counter 
do k=1,256 
Migweienu( in) jo xX Inieik, 7 )=243 
enddo 
enddo 
enddo 


Gems —1 256 
Go = 2536 
eee (epi i724 3 yx int (177) )=0 
enddo 
enddo 
return 
end 


OW 


* + + * 
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APPENDIX C. HEURISTIC SEARCH SOURCE CODE 
HS.FOR 


program hs2 

byte b image (256, 2 5c) 

integer sample image(256,256),clust image(256,250am 
counter/0/,s_peak_row(256),s_ peak col(25Gu— 
c peak row(256),c peak col(256),m,n, 
new_mark_row,new_mark col,min_m,min_n, 
minimam c n,Nc,row Gist,coleuase 

real*4 d(256,256),total cost(256),pnesene coc 
min,factor,hali dist, d elncee acre 
d_ clust total,Dc(256),minimum,row_dist_sq, 
col dist_sq 


open(unit=1,name='testper5.dat',status=‘old’, 
access='direct',recordsize=64) 


72* BEGIN Malt Peoguan~* ~ 
Read in test image. 
GiOr tai 2 a6 
read( 11 )(b Amage(ay )) 4-17 2 oo) 
enddo 


Convert input 256x256 byte image into 
a 256x256 integer image. 
call convert(b image,sample_ image,clust_image) 


Locate peaks in input image. 
do i=1,256 
ado. j— 1,256 
if(sample image(1,}]).NE.0)then 
COUNLEeEr-Counger+! 
S peak row(counter)=1 
c_peak_ row(counter)=1 
Ss peak _col(counter)=j 
c peak _col(counter)=j 
endif 
enddo 
enddo 


Calculate distances between peaks. 
call dist(counter,c peak row, ci peakuco la) 


7 0 


Seamer Manet wear am LOOp. 
do i=1,counter 
jJ=(counter+1)-1 
Cal INGererola( li, 3),4,COUuReEer ,min,s peak row, 
Spec CO lcm cake GOw,cC pecans coOl,Dc) 
call cost(}),Dc,present cost) 
EOrGMmMEOSE( J) =DEeCSENt Cost 
Daal: ), COLSIMCOST (4 ) 
feminat (ix, COealencOse( 321," )=', £8.4) 
enddo 


Gcumieman COst(tocal cost,counter) 
end 
end Main Pe Eogramn* ** 


RRR ION KR R Eee Ree er ee KKK KEKE KKK KKK Kee Kee Keke keke KK Kk 
Supbmpeutine CoOnvert(b image,sample image,clust_ image) 
byte Doimege (256, 256) 

DiPeeQck caloue wimage(256,2565),clust image(256, 256) 
real*4 factor 


faG@ror—243 ,.0/2 20 
do 1=1,256 
GOouj—1,256 
ee oeinacdc (iva). ll.0)then 
Sample image(i, 3)=b_image(i,})+256 
clust. image(i,}))=sample image(i,j) 
elseif(b_image(i,j).GE.0)then 
sample image(1i,j})=b image(1i,j) 
clust_image(i,j})=sample_ image(1i,)) 
endif 
sample image(i,j)=jnint(sample image(1i,j)/factor) 
clust_image(i,j)=jnint(clust_image(i,j)/factor) 
enddo 
enddo 
ae eet 
end 
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SUDTOUEINe GiSst(counter,c peak row,c peak _col,d) 
integer counter,c peak row(256),c peak _col(256) 
im ednan 4 enc 5 6.25.0) 


do m=1,counter 
do n=1,counter 
ec Nem) then 
ow Pe od@ernlode((c peak Trow(my—-c peak row(m))**2) 
Telodualempeak cOl(n)—c peak col(m))**2) ) 
endif 
enddo 
enddo 


ek 


CeEEUrRA 
end 
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subroutine centroid(i,j},d,counter,min,s peak row, 


. Ss peak_col,c peak row,¢ peak colypey 
integer 1,},counter,s peak row(256),s_ peak col( 25698 
* c peak row(256),c_peak col(256) 


real*4 d(256,256),m2in, De (25c) 


1f(37,EO counter) enen 
Det 7 =. 0 
else 
call min_dist(d,counter, min, ming ie 


Determine midpoint between two 

closest peaks. 

call new_marker(counter,min,min_m,min_n, 
* c peak row,c peak _col,hali diusee 
new_mark_row,new_mark_col) 


call clust_array ad}(counter,c peak iio, 
x c peak _col,new_mark_row, 
ss new _mark col,min_m,min_n) 


call dist_clust(j,counter,s peak row,s peak col 
‘ c peak row,c peak col,Dc) 


Calculate distances in adjusted imace. 
call dist (counter,c peak vrow,c peakmeo wig 
endif 
return 
end 


keke kee KR KKK KKK KKK KKK KK KKK KKK KKK KKK KKK KKK KKK KRKK 
subroutine min dist(d,counter,min,min mM, 02a 
integer counter,min_m,min_n 

real*4 <d¢256,256) mn 


do m=1,counter 
ao n=), countess 
if(d(m,n).NE.0)then 
min=d(m,n) 
Goer 
endif 
enddo 
enddo 
continue 


he 


Goem— i COUMECE 
do n=l,counter 
it (Genie 0.) then 
renee lame mL Tl ) teen 
min=d(m,n) 
min _m=m 
min n=n 
endif 
endif 
enddo 
enddo 
return 
end 
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subroutine new marker(counter,min,Min_m,min_n, 

2 GmpcomaloOwrc peak col, halt dist, 
x new _mark row,new_mark_col) 
Mimeger COuUNTeH, Mime m,Min n,c peak row(256), 

% Cupecanmeol( 7 56) ll, 1, New Mark row, néw mark col 
real*4 min,half dist 


m=min_m 
n=min_n 
min m=n 
min n=m 
Malt dist=—0.5*min 
if(c peak row(min_m).EQ.c peak row(min_n))then 
new mark row=c_peak row(min_m) 
New Meumecol—jnint( half distt 
Eloatvespeak col(min m))) 
Slice lt ecupecdianeo Mnlnem i. 2Osenpeak col(min n))then 
new mark» col=-c peak col(min_m) 
new_mark_row=jnint(half dist+ 
Eloat(c peak row(min m) )) 
Soom capeakeecol(minem).b..¢ peak cCol(min.n))then 
new_mark_row=jnint(float(c_peak_row(min_m) ) 
+ float(c peak row(min_n) 
-c peak _row(min_m))/2.0) 
lew Mogkecol—qnant (half wdist+ 
ioatCGepecakncoOl(mln im) ) ) 
eco ecm cc am ecOlm(Mimen)-Gl.c peak col(min n))then 
new_mark_row=jnint(float(c_peak_row(min_m) ) 
feulede(copean, COW( Min 1) 
-c peak row(min_m))/2.0) 
MeWolamiecOl—jnintk(t loat(c. peak —col(min_m) ) 
- -half dist) 
encour 
ime tard 
end 


* 


* 


* 


* 


US 
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subroutine clust_array_adj(counter,c peak row, 
x c peak _col,new_mark row, 
* new mark col,min m,™mingae 


integer counter,c_ peak _ row(256)/¢ peak com Zou, 
new_mark_row,new_mark col,min_m,min_n, 
A,B,C,D,E,F 


do n=1,couneer 
A=c_peak_row(n) 
B=c peak _row(min_m) 
C=c_ peak row(min_n) 
D=c_peak_col(n) 
E=c peak col(min im) 
F=c peak -col (mine 


if(((A.EQ.B).OR.(A.EO.C)).AND.((D.EQ.H) OR .( Deere eee 
c peak_row(n)=new_mark_row 
c peak col(n)=new_mark_col 
enait 
enddo 
perurn 
end 


Cc KREKEKREKKEEKKEKEKEKEKKEKEKREKEKKEKKEKKEEKKEKKKEKKKEEKX SKS 


subroutine dist clust(j,counter,sS peak row,s peak jeam 
: c peak row,c peace pe) 


integer j,counter,s peak row(256),Ss peakecoliZ5cy 
- C péak row(256)7C peak col( Zao 
row sdi1 St; cot sda se 
real*4 d_clust(256),d clust cotal) Der eiar 
: roW Gist sq7eol icteense 


d-clust fetal 00 


do n=1,counter 

row _dist=s_ peak_row(n)-c_ peak row(n) 

col dist=s_peak colin)-c  peakieoiGn 
row_dist_ sq=float(row_ dist*row dist) 

col dist sq=float(col dist*col dist) 
d_clust(n)=sqrt(row_dist_sq+col_dist_sq) 
d_ clust total=d clust(m)+G 7c is eewo rae 

enddo 


De(j])-Guclustleoea 


return 
end 
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Subroutine cost(j,Dc,present_cost) 

integer j,Nc 

real*4 Dce(256),present _ cost,Kn/1.0/,Kd/0.1/ 
Nc=j 

present _ cost=Kn*float(Nc)+Kd*Dc(j) 

ee cilltaml 

end 
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SuUDEOUGENIeC Mime ost( LOtal cost,counter ) 
integer counter,minimum_cn 
real*4 total _cost(256),minimum 
minimum=total cost(counter) 
do n=l1,counter 
mii bOtalm COSt(n) ._LE.minimum) then 
mMinimum=total cost(n) 
minimum _c_n=n 
endif 
enddo 
ete) eT) 
end 
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APPENDIX D: MULTIPLY.FOR SOURCE CODE 


MULTIPLY.FOR 

program multiply 

This program takes tracks.b, the fixed length 128x512 
byte output file from GIT, and converts 70 a aeem- 
256x256 byte array. Then the input LOFARGRAM, also a 
256x256 byte array, is multiplied with tracks.b. ie 
resultant image 1s identical to tracks, Deweeiee 
exception that each non-zero pixel location holds the 
value of beam00.bin at that i,j location. The ougeue 
file, OUTIMG.DAT, 1s used as input for the hough 
transform program. 


OA OY Oe) A557" 


byte bl image(128,512),b2_imege(2s0, 2508 
: b3 image(256,256),b4 image(256, 256) 


open(unit=1,name='tracks.b',status='old',access='direct ; 
*recordsize=128,maxrec=128) 
Oopen(unit=2,name='beam00.bin',status='old', 
*access='direct',recordsize=64,maxrec=256) 
open(unit=3,file='outimg.dat',status='new', 
*access='direct',recordsize=64,maxrec=256) 


do 1=1,128 
read( 1 1)(bl_ image. iy jj. )-1 oe) 
do (7 =145256 
b2 image((i*2)-1,))=bl_image(i,j) 
enddo 


dO] =75),, 572 
b2_ image(i*2,)-256)=bl_image(i, j) 
enddo 
enddo 


do 1=1,256 


read(2'1)(b3_i1mage(1)9)y ae 
enddo 


eS 


Gorr=1,256 
Gowns —15 256 
if(b2_ image(i,j).NE.0) then 
b4 image(i,j})=b3_image(i,}) 
else 
b4 image(i,j})=0 
endif 
enddo 
enddo 


do.i=1,256 
Writes seo 4 Tinage( 1, 3), }]=1,256) 
end do 
close(unit=1) 
Glose (Unit =Z ) 
close(unit=3) 


end 


ee 


APPENDIX E: HOUGH TRANSFORM (GRAY SCALE) SOURCE CODE 
HOUGH_G.FOR 


program hough g 


byte b _image(256),out_image(256,256),r image(Zacm 
integer counter,u(256),max,xgoneG 6 ue 

* gray level,h théeta( 256) herne 2eue 

* 1 _image( 256,256), aceum( 256, 2598 

* accum_norm(256,256),1t,1r,teSst)image@ s0 as 
real x0,y0,rho(256) ,x(256, 256) , p1/ 3a 

- theta(256),delta theta,rno mMaxyamemee 

* qélta eho; factor 


open(unit=1,name='outimg.dat',status='old', 


* access='direct',recordsize=64) 
open(unit=2,name='gttouthl6.dat',status='new', 

* access='direct',recordsize=64) 
open(unit=3,name='gttoutrl6.dat',status='new', 

. access='direct',recordsize=64) 
open(unit=4,name='gtthsl6.dat’,status='new', 

: access='direct',recordsize=64 ) 


***Begin Main Program*** 
Read in) bese image daca. 
doa y=17,256 
read(l'j) b_image 
dou = 1256 
1 image(i,j)=b_image(i) 
enddo 
enddo 


do 7=1 7256 
do-1=1 256 
if(i_image(i,j).LT.0) i image(i,})=1_image( 174) 
enddo 

enddo 
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Sar nougmh, clans tormm* ** 
Crmieinagedecum,cnieta,delta rho,rho 0,x0,y0) 
max=accum(1l,1) 
min=accum(1,1) 
do j=1,256 
Gomer 2516 
PCat |) wGlatlanx) mMax—accum( 1,3) 
Pacem ij). f bwin) mMin=accum( i, 7) 
enddo 
enddo 


Normalize accumulator array for display 
on IM-4000 Image Manager. 
factor=243.0/(max-min) 
Glog) = b, 2 510 
do i=1,256 
Qeaunmmemm( a, yi jl inlet aeoumi, 7) —-min)*factor ) 
if(accum_norm(i,j}]).GT.127) then 
b image(i)=accum_norm(i,]j)-256 
else 
Dmeragde,))—-accum NOrm( i, j) 
endif 
enddo 
write(2'j) b_image 
enddo 


Za—reecOnstructe from hough transftorm*** 
Tien (cCCulmnaCeCUMmmNOnN, Max, i.e theta,h rho,rho, 


* delea whe, rio. ,;x0,y0,x Int) 
eo 3=1,256 
do 1=1,256 


ae emer Gia 7) xX ine Ci,))=x int(i,j))-256 
r image(i)=x_int(i,j) 
enddo 
write(3'j) r_image 
enddo 
end 
oe eect he Program * ** 


kaekkk kkk kkk kkk kee K KKK KKK KKK KKK KKK KKK Ka KK KKK KKK KKK KKK Kk 
Subroutine h(i_image,accum,theta,delta_rho,rho_0,x0,y0) 
DiGegdoumilmIlmage(256,2596) ,accum(2560,2596),gray level 

real Peewee oo) <0, yO,pi1/3.1415926/ 


Mecmwenemeeneta trom 0 FO pi radians. 

delta theta=pi/float(256) 

@e-1=1,256 
theta(i)=float(i-1l)*delta theta 

enddo 


ws 


Q 


* 


* 
* 


* 


Fix the center of image as sehevon Lome 
rho _max=sqrt(float(256*256)+i toat( 256" 7am 
rho_O0=rho_max/2.0 
delta rho=rho max/filoaec( 256) 
xO=float(256/2) 
VO0=iloOae( 2567 2) 
de iy=13256 
Vey 
dO ixely 256 
gray level=1_ image(ix,iy) 
1f(gray tevel .bE.0 ) Ncomten.> 
X=1X 
do it=1,256 
Hough Transformation Equation. 
rho=(x-x0)*cos(theta(it))+(y0-y)*sin(theta(it) ) 
r=(rho+rho_0)/delta_rho 
1r= none) 
accum( it, 1r)=accum@t, wa ama eet e ve 


enddo 
continue 
enddo 
enddo 
return 
end 


WRAEKKKE RE KRKKKEKEKKREKRKEK KEKE KKEKR KKK KKH KKK RRR KK KS KK 


subroutine r(accum,accum_norm,max,h theta, horney aimee 
delta rho,rhoe.0,x0,y0, xine 

byte out_image(256, 256) 

integer accum(256,256),accum norm( 256,256), max, masa 
h_theta(256),h_rho(256),x_i1nt( 2567255) ue 
counter,test_image( 256,256) 

real rho(256),delta rho,rho 907 x0, 70) 402560) 2- 3. 
pi/3.1415926/ 


max th=jnint(0.16*max) 
Ceountenr—0 
ado j=15.2 20 
dou = 12516 
1f(accum(i,7).GE-maxcen ere 
GOUNTEr=COuneer. 
test_image(1,j )-accumeneni eu age 
Determine the theta and rho locations 
in the accumulator array with the 
highest count. vatue- 
h theta(counter)=i 
h_ rho(counter)=j 
endif 
enddo 
enddo 
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do 1=1,256 
Ger y=1 256 
if(test_image(i1,j]).GT.127)then 
Quem vac owivg))—teStuiipage{ 1, ])—=296 
else 
out image(i,j)=test_image(i,}j) 
end if 
enddo 
enddo 


gio i= 256 
write(4'1)(out_image(i,j),j=1,256) 
enddo 


do j=l1,counter 
Riot} )—tmeuner |) "aelta rho-rho 0 
enddo 


Ge j=), counter 
Geoe1=1,256 
RECOllcoeSUchren EQuation. 
rae, ie Oe 
Conewmpreeww- 1) sin(( ho ehetat ])-1)*pi/256) ) / 
COc(iaeeneta( )})—1)*p1/256) 
x_int(i,j)=jnint(x(i,j)) 
enddo 
enddo 


SOnj-1, 256 
do l=1,counter 
Hele) — eee} |) 
enddo 


do m=l1,counter 
dion k=17.25.6 
Paes Oou (moe 1 nitiikyes))=243 
enddo 
enddo 
enddo 


Gor y= 1,256 
do 1=1,256 
Memes, J) NE. 243) xX int(i,;))=0 
enddo 
enddo 
eee al 
end 
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